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The  steady-state  solution  of  the  system  of  equations  consisting  of  the  full  Navier-Stokes  equa¬ 
tions  and  two  turbulence  equations  has  been  obtained  using  a  multigrid  strategy  on  unstruc¬ 
tured  meshes.  The  flow  equations  and  turbulence  equations  are  solved  in  a  loosely  coupled 
manner.  The  flow  equations  are  advanced  in  time  using  a  multi-stage  Runge-Kutta  time  step¬ 
ping  scheme  with  a  stability  bound  local  time-step,  while  the  turbulence  equations  are 
advanced  in  a  point-implicit  scheme  with  a  time-step  which  guarantees  stability  and  positivity. 
Low  Reynolds  number  modifications  to  the  original  two-equation  model  are  incorporated  in  a 
manner  which  results  in  well  behaved  equations  for  arbitrarily  small  wall  distances.  A  variety 
of  aerodynamic  flows  are  solved  for,  initializing  all  quantities  with  uniform  freestream  values. 
Rapid  and  uniform  convergence  rates  for  the  flow  and  turbulence  equations  are  observed. 
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1.  INTRODUCTION 


The  use  of  unstructured  meshes  has  become  more  widespread  in  recent  years  due  to  the 
ease  with  which  complex  geometries  can  be  handled  and  the  possibility  of  enhancing  the  solu¬ 
tion  accuracy  and  efficiency  through  adaptive  meshing  techniques.  To  date,  most  of  the 
successes  of  unstructured  mesh  techniques  have  been  in  computing  inviscid  flows  in  two  and 
three  dimensions  over  arbitrary  geometries.  However,  more  recently,  solutions  of  the  Navier- 
Stokes  equations  on  unstructured  meshes  have  been  reported  [1,2,3, 4,5].  The  main  obstacles  to 
efficiently  computing  high-Reynolds-number  flows  on  unstructured  meshes  are  due  to  the 
required  grid  stretching  and  the  turbulence  model.  For  high-Reynolds-number  flows  over 
streamlined  bodies,  viscous  effects  are  confined  to  thin  boundary-layer  and  wake  regions, 
which  can  only  be  resolved  efficiently  using  high  aspect  ratio  elements.  One  approach  [3,5]  is 
to  fit  a  thin  local  mesh  of  structured  high  aspect  ratio  quadrilaterals  in  the  viscous  regions,  and 
fill  the  remainder  of  the  domain  with  an  unstructured  mesh.  The  other  approach  consists  of 
filling  the  entire  domain  with  an  unstructured  mesh  which  contains  highly  stretched  triangular 
elements  in  the  viscous  regions  [4].  In  this  work,  the  latter  approach  has  been  pursued,  in  the 
interest  of  developing  a  more  general  method  capable  of  dealing  with  a  wider  variety  of  flows, 
such  as  flows  with  confluent  boundary  layers,  or  mixing  wakes,  and  also  to  enable  the 
straight-forward  implementation  of  adaptive  meshing  techniques  throughout  all  regions  of  the 
flow-field.  The  numerical  scheme  must  therefore  be  formulated  such  that  the  accuracy  and 
convergence  are  not  seriously  affected  by  the  presence  of  highly  stretched  triangular  elements. 

The  most  commonly  employed  turbulence  models  for  compressible  flow  calculations  are 
of  the  algebraic  mixing-length  type  [6].  These  models  have  been  shown  to  produce  good 
results  for  attached  turbulent  boundary  layers  and  mildly  separated  flows  using  structured 
meshes,  and  have  also  been  implemented  for  non-trivial  geometries  on  unstructured  meshes  [7]. 
Although  such  models  can  be  made  inexpensive  and  computationally  robust  even  in  the  context 
of  unstructured  meshes,  they  lack  the  generality  required  for  dealing  with  completely  arbitrary 
geometries,  and  their  ability  in  predicting  flows  with  multiple  confluent  shear  layers  and  large 
amounts  of  separation  is  at  best  limited.  Two  equation  models,  on  the  other  hand,  offer  the 
possibility  of  dealing  with  the  more  complicated  flows  which  are  often  associated  with  the 
complex  geometries  for  which  unstructured  meshes  are  so  well  suited.  In  principle,  the  imple¬ 
mentation  of  such  models  on  unstructured  meshes  can  be  accomplished  in  a  straight-forward 
fashion,  simply  by  discretizing  and  integrating  the  turbulence  equations  in  a  manner  analogous 
to  that  employed  for  the  mean  flow  equations.  However,  field-equation  turbulence  models  have 
often  proved  to  be  extremely  difficult  to  integrate  to  steady-state,  exhibiting  stiff  or  unstable 
numerical  behavior  in  regions  very  close  to  the  wall,  as  well  as  in  the  far-field.  The  use  of 
multigrid  to  solve  the  turbulence  equations  has  recently  been  reported  by  several  authors  [8,9], 
using  a  Ni-type  scheme  on  structured  meshes.  In  this  work,  a  multigrid  strategy  which  has 
previously  been  developed  for  the  Euler  and  Navier-Stokes  equations  on  unstructured  meshes 
[4,10]  is  extended  to  solve  for  the  two  turbulence  equations  as  well. 

2.  GOVERNING  EQUATIONS 

The  governing  equations  are  obtained  by  Favre  averaging  the  Navier-Stokes  equations, 
and  modeling  the  Reynolds  stress  and  heat  flux  terms  by  the  Boussinesq  assumption.  In  con¬ 
servative  form,  these  equations  are  written  as 

dw  dfc  fa  d/v 
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where  w  is  the  solution  vector  and  fc  and  gc  are  the  cartesian  components  of  the  convective 
fluxes 
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In  the  above  equations,  p  represents  the  fluid  density,  u  and  v  the  x  and  y  components  of  fluid 
velocity,  E  the  total  energy,  and  p  is  the  pressure  which  can  be  calculated  from  the  equation  of 
state  of  a  perfect  gas 

p  =  (Y-l)p  E  -  (3) 

The  viscous  fluxes  /„  and  gv  are  given  by 

fo  1  fo  1 


uoa+va„-q. 


where  <r  represents  the  stress  tensor,  and  q  the  heat  flux  vector,  which  are  given  by 
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P  represents  the  molecular  viscosity,  and  p,  denotes  the  turbulent  eddy  viscosity,  which  must 
be  computed  by  a  suitable  turbulence  model.  Pr  is  the  laminar  Prandtl  number,  which  is  taken 
as  0.7  for  air,  Pr,  is  the  turbulent  Prandtl  number,  taken  as  0.9,  and  7  is  the  ratio  of  specific 
heats  of  the  fluid. 

The  high  Reynolds  number  k-e  turbulence  model  originally  described  by  Launder  and 
Spalding  [11],  can  similarly  be  written  as 

Bw  ,  d/c  ,  dgc  d/v  fa  , 

T  +  IT  + 17  =  "aT  + 17  +  *  <10> 

where  w ,  fe  and  gc  are  now  given  by 


The  diffusive  fluxes  fv  and  gv  are  given  by 
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and  the  source  term  h  is  given  by 
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where  the  production  term  P  and  the  term  S  in  two  dimensions  are  given  by 

p  =  j(u?  +  Vy2  -  u,vy)  +  (Uy  +  V,)2  (14) 

S  =  u,  +  v. 

The  eddy  viscosity  is  calculated  as 

pCJi2 

*  =  e  (15> 

and  k  also  appears  in  the  normal  stresses  in  equation  (S).  The  constants  appearing  in  the 
above  equations  are  given  the  standard  values  recommended  in  [11],  i.e. 

CV  =  0.09  or*  =  1.0  ae  =  1.3  Cx  =  1.44  C2  =  1.92  (16) 

These  equations  are  coupled  to  the  governing  equations  for  the  mean  flow  and  exhibit  a  similar 
structure.  Therefore,  a  single  system  of  equations  which  simultaneously  governs  the  flow  and 
turbulence  quantities  may  be  written  as 

+  +  A  /i  tv 
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dt  dx  dy  dx  dy 

where  the  solution  vector  and  the  source  term  are  now  given  by 


[pe  J  M.'’  -  jSpt  - (* 

and  the  flux  definitions  follow  from  equations  (2),(4),(1 1),  and  (12). 

The  solution  procedure  consists  of  discretizing  these  equations  in  space  on  an  unstruc¬ 
tured  mesh,  and  then  integrating  the  discretized  equations  in  time  until  the  steady-state  solution 
is  obtained.  The  basic  strategy  pursued  in  this  work  involves  the  use  of  a  finite-element  Galer- 
kin  discretization  technique,  in  conjunction  with  an  unstructured  multigrid  integration  technique 
to  solve  for  the  steady-state.  Although  all  six  equations  of  the  governing  system  are  solved 
simultaneously  in  the  multigrid  strategy,  the  flow  equations  are  only  loosely  coupled  to  the  tur¬ 
bulence  equations  (through  the  value  of  p,),  and  we  choose  to  employ  somewhat  different  base 
grid  solvers  for  the  flow  equations  and  the  turbulence  equations. 

3.  SPATIAL  DISCRETIZATION 


The  equations  governing  the  mean  flow  are  discretized  using  a  Galerkin  finite-element 
approach  [4].  The  flow  variables  are  stored  at  the  vertices  of  the  triangles.  The  convective 
fluxes  are  computed  at  the  vertices  of  the  triangles  and  assumed  to  vary  linearly  over  the  tri¬ 
angular  elements.  For  the  viscous  terms,  the  flow  variables  themselves  are  assumed  to  vary 
linearly  over  the  triangular  elements  of  the  mesh,  and  the  required  velocity  gradients  in  the 
expression  for  the  viscous  stresses  are  thus  computed  at  the  centers  of  the  triangular  elements. 
Additional  artificial  dissipation  terms  are  required  to  ensure  the  stability  of  the  convective 
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terms  and  these  are  constructed  as  a  blend  of  a  Laplacian  and  bihaimonic  operators  in  the  con¬ 
served  variables,  designed  to  ensure  second-order  accuracy  throughout  the  flow-field,  except  in 
the  vicinity  of  a  shock  where  first-order  accuracy  is  recovered.  For  the  turbulence  equations, 
the  diffusive  terms  are  similarly  discretized  using  a  Galerkin  finite-element  approach,  assuming 
linear  variations  of  the  conserved  variables  over  the  triangular  elements.  The  velocity  gradients 
in  the  source  terms  are  also  constructed  assuming  linear  elements.  The  convective  terms,  how¬ 
ever,  are  constructed  using  first-order  upwinding.  Although  only  first-order  accurate,  this 
approach  is  employed  since  it  helps  ensure  stability  and  positivity  of  the  conserved  variables 
throughout  the  integration  procedure,  as  will  be  shown.  Furthermore,  in  regions  where  convec¬ 
tion  is  small  compared  to  the  diffusion  terms  or  the  source  terms,  such  as  in  the  logarithmic 
law  of  the  wall  region,  the  scheme  reverts  to  second  order  accuracy.  In  future  work  however,  a 
second-order  accurate  implementation  of  the  convective  terms  may  be  pursued. 


4.  INTEGRATION  SCHEME 


The  discretized  mean  flow  equations  are  integrated  in  time  using  an  explicit  five-stage 
Runge-Kutta  time-stepping  scheme,  where  the  convective  terms  are  evaluated  at  every  stage, 
and  the  dissipative  terms  are  only  evaluated  at  the  first,  third  and  fifth  stages.  This  scheme, 
which  has  previously  been  described  [4,12],  has  been  particularly  devised  to  ensure  rapid 
damping  of  high-frequency  errors,  and  is  thus  well  suited  to  drive  the  multigrid  algorithm. 
Convergence  is  accelerated  by  the  use  of  local  time-stepping,  and  implicit  residual  averaging. 
In  principle,  the  turbulence  equations  may  be  integrated  in  time  using  the  same  explicit 
scheme.  However,  the  presence  of  source  terms  imposes  a  further  time-step  restriction.  If  the 
flow  equations  and  turbulence  equations  are  integrated  in  a  fully  coupled  manner,  the  minimum 
local  time-step  from  the  flow  and  turbulence  equations  must  be  employed.  In  regions  where  the 
source  terms  dominate,  this  may  lead  to  slow  convergence.  If,  on  the  other  hand,  the  flow 
equations  and  turbulence  equations  are  integrated  in  an  uncoupled  explicit  manner,  the  tur¬ 
bulence  equations  may  significantly  lag  the  flow  equations  and  thus  inhibit  convergence  to  the 
steady-state  solution.  In  order  to  advance  the  turbulence  quantities  at  the  same  rate  as  the  flow 
equations,  the  source  terms  must  be  treated  implicitly.  However,  rather  than  simply  treat  the 
source  terms  implicitly,  the  system  of  turbulence  equations  is  integrated  in  a  point-implicit 
manner.  Thus  we  rewrite  the  discretized  turbulence  equations  as 


aw. 

=  l?(w.)  +  H(Wi)  (19) 

At 

where  R  (w, )  represents  the  discretized  convective  and  diffusive  terms,  which  depend  on  the 
values  of  w  at  i  and  at  neighboring  nodes,  and  //(*♦',•)  represents  the  discretized  source  terms, 
which  only  depend  on  the  values  of  w  at  i.  The  above  equation  is  then  linearized  about  the 
values  at  i  which,  upon  solving  for  Aw,-  yields 


Ah';  = 


(20) 


The  Runge-Kutta  scheme  described  above  is  now  replaced  by  a  multi-stage  implicit  scheme, 
where  the  qth  stage  is  given  by 


w 


<«>  =  W«J> 


a,  A/ 


dR 

dw 


dH 

dw 


(flfr-'V.)  +  //(,"’V.)j 


(21) 


where  the  a,  denote  the  Runge-Kutta  coefficients  for  the  qth  stage,  and  Ar  is  the  local  time 
step.  In  this  manner,  the  high-frequency  damping  characteristics  of  the  original  scheme  are 
approximated,  while  the  time-step  restriction  due  to  the  source  terms  is  alleviated.  The  precise 


-5- 


value  of  the  local  time-step  At  employed  is  one  which  guarantees  stability  as  well  as  positivity 
of  the  turbulence  quantities. 

5.  STABILITY  AND  POSITIVITY  CONSIDERATIONS 

One  method  to  guarantee  stability  of  the  system  is  to  ensure  that  the  matrix  to  be  inverted 
is  diagonally  dominant.  This  is  not  a  necessary  condition  for  stability,  although  it  is  sufficient. 
This  can  obviously  be  achieved  by  choosing  At  to  be  sufficiently  small.  However,  the  reason 
for  employing  a  point-implicit  approach  now  becomes  apparent.  Since  the  two  turbulence  equa- 

dR 

tions  are  only  coupled  through  their  source  terms,  the  —  matrix  is  diagonal.  The  contribution 

from  the  diffusive  terms  is  strictly  negative,  as  well  as  that  from  the  first-order  upwinded  con¬ 
vective  terms.  Hence,  these  terms,  when  subtracted  from  the  diagonal  of  the  matrix  to  be 
inverted,  increase  the  diagonal  dominance,  and  hence  permit  the  use  of  a  larger  time-step.  The 
maximum  value  of  At  is  found  by  equating  each  diagonal  element  to  its  corresponding  off- 
diagonal  element  in  the  coefficient  matrix.  The  actual  value  employed  for  the  time-step  is 
taken  as  the  minimum  between  the  two  values  obtained  by  the  diagonal  dominance  test,  and 
the  value  determined  by  local  stability  analysis  for  an  explicit  scheme  in  the  absence  of  source 
terms. 

Physically,  k  and  e  represent  quantities  which  must  remain  non-negative.  Thus  a  further 
time-step  restriction  is  required  to  ensure  positivity.  For  a  simple  2x2  system,  this  can  easily  be 
derived  analytically.  Thus,  we  require  that  the  new  update  to  the  turbulence  variables  be  such 
that 

w  +  Aw  >  aw  0<a<l 

or,  when  Aw  <  0 

lAw  I  <  (1-a)  w  0  <  a  <  1  (22) 

Substituting  into  equation  (20),  and  using  Cramer’s  rule  to  evaluate  the  inverse  of  the  2x2 
matrix,  we  obtain  two  quadratic  inequalities  for  At,  i.e.  one  for  positivity  of  jfc,  and  one  for  e. 
The  time  step  is  then  limited  by  the  smallest  positive  root  of  the  two  quadratic  equations. 

6.  MULTIGRID  STRATEGY  AND  STEADY-STATE  CONSIDERATIONS  FOR  THE 
k  -  e  EQUATIONS 

A  multigrid  strategy  is  employed  to  accelerate  the  solution  of  the  system  of  mean  flow 
and  turbulence  equations  to  steady-state.  In  the  context  of  unstructured  meshes,  multigrid  may 
be  applied  by  generating  a  sequence  of  non-nested  coarse  and  fine  meshes,  and  transferring  the 
variables,  residuals  and  corrections  back  and  forth  between  the  various  meshes  using  linear 
interpolation.  The  patterns  for  interpolating  between  non-nested  unstructured  meshes  are  deter¬ 
mined  in  a  preprocessing  stage,  using  an  efficient  search  algorithm.  The  present  multigrid  stra¬ 
tegy  has  previously  been  described  in  detail  for  the  Euler  and  Navier-Stokes  equations  [4,10], 
and  thus  will  not  be  repeated  here.  In  previous  multigrid  applications  for  turbulent  flows  using 
algebraic  models  on  structured  and  unstructured  meshes  [7,13],  the  eddy  viscosities  are  only 
computed  on  the  finest  grid,  and  interpolated  to  the  coarser  meshes.  Since  the  eddy  viscosity 
represents  the  main  coupling  between  the  flow  equations  and  the  turbulence  equations,  a  simi¬ 
lar  approach  has  been  adopted  in  the  present  context,  thus  ensuring  a  more  accurate  representa¬ 
tion  of  this  quantity  on  all  grid  levels.  However,  since  the  eddy  viscosity  is  only  computed  on 
the  finest  grid,  it  is  effectively  held  constant  throughout  an  entire  multigrid  cycle,  and  the 
source  terms  must  be  linearized  accordingly.  Making  use  of  equations  (11)  and  (13),  the  linear¬ 
ization  of  the  source  terms  on  all  grids  is  therefore  taken  as 
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At  this  point  it  is  worth  commenting  on  what  types  of  errors  may  be  expected  to  be  handled 
efficiently  by  a  multigrid  strategy.  Multigrid  is  an  effective  device  for  relieving  the  spatial 
stiffness  associated  with  a  set  of  discretized  equations,  which  is  achieved  by  time-stepping  on 
coarser  grids.  The  turbulence  equations  contain  spatial  terms  such  as  convection  and  diffusion, 
but  the  source  terms  are  purely  local  terms.  In  fact,  in  the  absence  of  convection  and  diffusion, 
the  equations  become  completely  uncoupled  in  space,  and  a  property  formulated  multigrid 
algorithm  should  yield  vanishingly  small  corrections  in  such  a  case.  Thus,  it  is  important  for 
the  base  grid  solver  to  efficiently  eliminate  errors  associated  with  these  terms.  From  another 
point  of  view,  if  a  purely  explicit  scheme  were  employed,  a  time-step  restriction  would  arise 
from  the  convection,  diffusion  and  source  terms.  While  the  first  two  restrictions  are  relaxed 
when  going  to  coarser  grids,  the  latter  remains  the  same  on  all  grid  levels,  effectively  prevent¬ 
ing  the  use  of  large  time-steps  on  coarse  grids  and  severiy  limiting  the  overall  rate  of  conver¬ 
gence.  The  use  of  a  point-implicit  scheme,  therefore,  relieves  any  such  restrictions,  and  results 
in  overall  convergence  rates  similar  to  that  achieved  with  the  mean  flow  equations. 


At  steady-state,  the  turbulence  equations  do  not  necessarily  exhibit  a  unique  solution.  In 
regions  where  the  production  term  \i,P  vanishes,  k  =  0,  e  *  0  is  an  obvious  solution  which  can 
be  found  by  inspection  of  equations  (10)  and  (13).  However,  the  eddy  viscosity,  which  is  given 
by  equation  (15),  becomes  a  ratio  of  two  vanishing  quantities,  and  is  thus  undefined.  The  time 
dependent  turbulence  equations  however,  are  not  ill-posed.  On  the  contrary,  the  value  of  the 
constant  C2  has  been  carefully  chosen  to  ensure  that  k,  e  and  p,  all  vanish  asymptotically  for 
an  isotropic  decaying  turbulence.  For  an  isotropic  turbulence,  all  spatial  terms  as  well  as  the 
production  term  vanish,  and  equations  (10)  and  (13)  reduce  to 


dk 

dt 


=  -  e 


(24) 


de. 

dt 


=  -  Co—- 


Solution  of  this  system  yields  the  following  asymptotic  behavior 


*  =  r 


i 

c2-i 


cri 


e  =  r  “ '  y  ~  t  r— (25) 

which,  for  the  current  value  of  1  <  C2  <  2  indicates  that  all  quantities  vanish  for  large  L 
Hence,  in  order  to  converge  to  the  appropriate  steady-state  solution,  it  is  important  for  any 
numerical  scheme  to  respect  the  relative  asymptotic  time  behavior  of  k  and  e  throughout  the 
convergence  process.  For  the  base  grid  solver,  this  is  achieved  by  employing  the  maximum 
time-step  for  the  system  of  two  turbulence  equations  which  ensures  stability  and  positivity;  let¬ 
ting  k  or  t  become  negative,  or  taking  too  large  time  steps  and  subsequently  limiting  the 
updated  values  of  k  or  e  invariably  leads  to  unrealistic  values  of  p,  in  the  far-field.  Within  the 
multigrid  strategy,  corrections  interpolated  back  from  the  coarser  grids  may  cause  k  or  e  to 
become  negative.  Rather  than  limit  these  corrections,  they  are  simply  omitted  at  any  point 
where  positivity  cannot  otherwise  be  guaranteed.  In  this  manner,  the  (point-wise)  time  con¬ 
sistency  is  not  violated,  and  the  overall  effect  is  simply  to  lag  such  points  by  the  effective 
coarse  grid  time  step.  An  alternate  approach  would  be  to  recompute  the  coarse  grid  corrections 
employing  a  smaller  time  step  which  guarantees  positivity.  However,  due  to  the  recursive 


c,-i 
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nature  of  multigrid,  this  represents  a  non-trivial  task. 


7.  BOUNDARY  AND  INITIAL  CONDITIONS 

The  derivation  of  the  above  turbulence  transport  equations  is  made  under  the  hypothesis 
of  a  large  Reynolds  number  flow.  Thus,  in  regions  close  to  the  wall,  such  as  in  the  viscous 
sublayer  where  molecular  effects  become  important,  these  equations  are  not  valid.  In  order  to 
avoid  integrating  the  turbulence  equations  in  these  region  we  make  use  of  wall  functions.  In 
this  approach,  the  governing  equations  for  the  flow  and  the  turbulence  are  integrated  up  to  a 
distance  y  =  yK  away  from  the  wall.  The  flow  in  the  remaining  region  0  <  y  <  yn  is  assumed  to 
obey  the  law  of  the  wall  i.e. 


U  1 .  pUj 
—  =  —In - 

Ux  K  ]i 


+  5.5 


(26) 


At  each  time-step  in  the  solution  procedure  of  the  governing  equations,  an  estimate  of  the  velo¬ 
city  U  at  y  =  yn  is  obtained.  From  this,  the  value  of  the  wall  shear  stress  can  be  obtained  by 
solving  equation  (26)  implicitly  for  Ux  (  using  a  Newton-Raphson  method).  This  estimate  of 
the  wall  shear  stress  is  then  employed  as  a  boundary  condition  on  the  momentum  equation  for 
the  mean  flow,  and  results  in  Dirichlet  wall  boundary  conditions  for  k  and  e.  In  practice  the 
point  y  =  y„  is  very  close  to  the  wall  so  that  it  may  be  approximately  placed  on  the  wall,  and 
the  boundary  conditions  at  y  =y„  may  be  imposed  at  the  wall  surface.  For  the  momentum 
equation,  this  results  in  a  wall  slip  velocity  U  =  U(yH). 


In  the  far-field,  k  and  e  are  assigned  freestream  values  at  inflow  boundaries,  and  simple 
extrapolation  is  employed  at  outflow  boundaries.  Initial  conditions  on  k  and  e  are  obtained  by 
imposing  a  level  of  freestream  turbulence  from  which  k  is  determined,  and  e  is  evaluated  from 
equation  (15)  in  order  to  produce  a  low  value  of  freestream  eddy  viscosity  (p,  <  1).  However, 
since  the  present  formulation  results  in  a  small  value  of  p,  in  all  regions  of  the  flow  field 
where  production  is  negligible,  the  converged  solution  is  relatively  insensitive  to  the  initial 
values  of  k  and  e.  The  mean  flow  equations  are  initialized  using  uniform  freestream  flow  con¬ 
ditions,  and  applying  the  tangential  slip  velocity  boundary  condition  (as  for  an  inviscid  flow). 
Throughout  the  integration  process,  the  wall  shear  stress  obtained  from  equation  (26),  which  is 
fed  back  into  the  momentum  equation,  retards  the  flow  near  the  wall,  thus  creating  a  boundary 
layer  profile. 


8.  RESULTS  USING  WALL  FUNCTIONS 

Two  attached  flow  cases  have  been  computed  using  the  multigrid  implementation  of  the 
high-Reynolds-numbcr  turbulence  model  described  above.  The  first  case  consists  of  transonic 
flow  past  an  RAE  2822  airfoil.  The  freestream  Mach  number  is  0.729,  the  incidence  is  2.31 
degrees,  and  the  Reynolds  number  is  6.5  million.  This  case  (case  6)  has  been  well  documented 
both  experimentally  [14]  and  computationally  [7,13,15]  on  structured  and  unstructured  meshes. 
The  mesh  employed  for  this  case  is  depicted  in  Figure  1.  It  contains  12,823  vertices  and  exhi¬ 
bits  a  normal  spacing  of  KT4  chords  at  the  airfoil  surface,  which  positions  the  first  point  off  the 
wall  in  the  logarithmic  law  of  the  wall  region.  The  computed  Mach  contours  for  this  case  are 
shown  in  Figure  2,  while  the  resulting  eddy  viscosity  distribution  is  given  in  Figure  3.  A 
smooth  distribution  of  eddy  viscosity  throughout  the  boundary-layer  and  wake  regions,  and 
vanishingly  small  values  in  the  inviscid  regions  of  flow  are  observed.  The  computed  surface 
pressure  distribution  is  compared  with  experimental  data  [14]  in  Figure  4,  showing  good 
overall  agreement.  The  convergence  rate  of  the  system  of  equations  is  depicted  in  Figure  5,  by 
plotting  the  RMS  average  of  the  density  residual,  and  the  residual  of  p k  throughout  the  flow 
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field,  versus  the  number  of  multigrid  cycles.  As  can  be  seen,  the  flow  equations  and  turbulence 
equations  converge  with  the  same  asymptotic  rates.  The  residuals  are  reduced  by  roughly  6 
orders  of  magnitude  over  200  cycles,  yielding  an  overall  convergence  rate  of  0.93 

The  second  case  involves  flow  over  a  high-lifting  four-element  airfoil.  This  case  is  useful 
in  demonstrating  the  complex  geometries  and  resulting  flow-fields  which  can  be  handled  by  the 
present  methodology.  The  mesh  employed  is  depicted  in  Figure  6.  It  contains  a  total  of  51,100 
vertices  and  a  normal  spacing  of  2x  KT*  chords  off  the  wall  for  each  airfoil  element.  The  com¬ 
puted  Mach  contours  are  shown  in  Figure  7,  while  the  resulting  eddy  viscosity  distribution  is 
given  in  Figure  8.  The  ease  with  which  multiple  wakes  and  confluent  boundary  layers  may  be 
handled  by  the  present  approach  is  evident  from  the  figures.  The  computed  surface  pressure 
distribution  is  seen  to  compare  favorably  with  experimental  wind-tunnel  data  from  [16]  in  Fig¬ 
ure  9.  It  should  however  be  pointed  out  that  such  favorable  agreement  is  in  large  part  due  to 
the  attached  nature  of  the  flow.  The  multigrid  convergence  rates  of  the  density  equation  and  the 
k  equation  are  depicted  in  Figure  10,  where  both  equations  are  seen  to  achieve  approximately 
the  same  asymptotic  rates,  decreasing  by  4  orders  of  magnitude  over  300  cycles. 

9.  LOW  REYNOLDS  NUMBER  TURBULENCE  MODEL  MODIFICATIONS 

While  the  use  of  the  high-Reynolds-number  turbulence  equations  in  conjunction  with  wall 
functions  is  useful  for  a  large  class  of  wall  bounded  flows,  it  is  nevertheless  limited  to  flows 
where  a  logarithmic  law  of  the  wall  region  exists,  and  is  thus  strictly  not  valid  for  separated 
flows.  An  alternative  approach  is  to  modify  the  turbulence  equations  in  order  to  account  for 
low-Reynolds-number  effects.  Many  such  modifications  have  been  proposed  over  the  years 
with  varying  degrees  of  success  [17].  One  common  feature  of  all  such  modifications  is  that 
they  have  proved  exceedingly  difficult  to  integrate  numerically  very  close  to  the  wall.  The  aim 
of  the  present  work  is  to  develop  an  efficient  and  robust  technique  for  integrating  such  models, 
rather  than  reformulating  or  advocating  any  one  model  in  particular.  With  this  in  mind,  we 
chose  to  implement  the  simplest  possible  low-Reynolds-number  model  that  has  been  demon¬ 
strated  to  produce  good  results  for  simple  problems,  with  possible  extensions  to  more  complex 
models  should  the  original  version  prove  inadequate  for  more  complicated  flows.  To  this  end, 
the  modifications  proposed  by  Speziale,  Abid  and  Anderson  [18]  have  been  implemented.  The 
modified  turbulence  equations,  now  given  in  vector  form,  can  be  written  as 


3pe 

3/ 


'  ' 

^  +  «.V(p*)  =  V.  (n+  —  )V*  +  -  ~S p*  -  pe 

dt  sk  3 

r  v 


+  n.V(pe)  =  V|(p+^)Vej  +  C^P  -  fSp*)J  -  C/zPy 

^  _  pcy  /n*2 

/»  =  (1  +  #)tanh(*r) 


/2  =  [  1  -  exp< -jj)  ]2 


with  boundary  conditions  at  the  wall  given  by 


k  =  0 
e  = 

P  dy2 

and  employing  the  following  values  for  the  constants 


(27) 


(28) 


a*  =  1.36 


ae  =  1.36 
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C ^  =  0.09 
where  Re, 


=  P 


cl  =  1.44  C2  =  1.83  [1  -  |exp(-^-)2  ] 


k2 

—  is  the  turbulence  Reynolds  number.  As  can  be  seen,  no  extra  source  terms 
\u 


are  introduced,  and  only  two  damping  functions  are  required.  This  model  is  similar  in  form  to 
the  Lam-Bremhorst  model  [19],  with  the  notable  difference  that  all  damping  functions  depend 
solely  on  y+.  The  evaluation  of  such  functions  requires  the  knowledge  of  the  distance  of  each 
point  from  the  closest  wall.  In  the  context  of  unstructured  meshes,  this  information  can  be  con¬ 
structed  through  the  use  of  a  generalized  distance  function,  as  outlined  by  Barth  [20].  As  with 
most  low-Reynolds-number  turbulence  models,  the  current  form  of  the  model  has  been 
reported  to  be  extremely  stiff  in  near-wall  regions,  generally  requiring  the  prescription  of  initial 
profiles  in  k  and  e  in  order  to  guarantee  convergence  to  steady-state.  Such  techniques  are  con¬ 
sidered  impractical  for  complex  aerodynamic  flows,  and  thus  a  more  robust  solution  strategy 
has  been  pursued.  The  difficulties  associated  with  the  near-wall  regions  can  be  assessed  by 
inspection  of  equations  (27).  When  the  wall  boundary  condition  it  =  0  is  substituted  into  the  e 
equation,  it  is  seen  to  result  in  a  singularity,  since  k  appears  in  the  denominator  of  this  equa¬ 
tion.  Since  /  2  also  vanishes  at  the  wall,  this  singularity  is  in  principle  removable.  However, 
the  numerical  integration  of  the  e  equation  in  its  present  form  will  only  be  well  behaved  if  /2 
and  k  have  the  same  asymptotic  behavior  near  the  wall  throughout  the  integration  procedure, 
thus  the  need  for  startup  profiles.  The  approach  taken  in  this  work  is  to  remove  the  singularity 
by  solving  for  a  new  variable  defined  as 


k=kf2  or  k  =  -~  (29) 

/  2 

Upon  substituting  this  expression  into  equation  (27),  and  using  the  chain  rule  to  evaluate  the 
gradient  operators,  one  obtains  the  new  set  of  equations 


^  +  «.V(pf)  =  V. 


(p  +  —  )V£ 

Sk 

1  7 

+  J-  [»,P  -jSpkfz  -  pe] 

+  7-  [  [  V(p  +  — )  -  pu  ]  .V/2  +  (p  +  —  )V2/2  ]  * 

J  2  Sk  *k 

+  7-  [  2  (n  +  )  V/2  .V*  ] 

J  2  Sk 


(30) 


3pe 

dt 


+  K.V(pe)  =  V. 


(p  +  —  )Ve 

•^e 


\i,p 


+  C,  7-7  -  f  C.Spe  -  C2pf 
J2  k  *  k 


While  the  k  equation  now  looks  rather  complicated,  the  new  terms  are  only  significant  in  the 
region  /2«  1,  and  in  fact,  although  all  terms  have  been  included,  only  the  £v2/2  term  has 
been  found  to  have  a  significant  effect  on  the  overall  solution.  At  the  wall,  we  have 


f  2  =  0,  V/2  =  0,  V2/2>0 

as  well  as 


p,  =  0  P  =  0  S  =  0  k  =  0 

The  boundary  condition  k  =  0  implies  H  bounded  at  the  wall.  Since  f2  which  appears  as  a 
denominator  in  the  right-hand  side  of  the  i  equation  vanishes,  we  require  the  non  vanishing 
terms  in  the  numerator  to  sum  to  zero,  thus  yielding  the  condition 


k  = 


pe  . 

pV2/2 


(31) 
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Upon  substituting  this  expression  into  the  e  equation,  with  u  =  0,  p,  =  0,  and  S  =  0  one  obtains 
a  modified  Helmoltz  equation  for  e  at  the  wall  in  the  steady-state 

v[pVe]  -  C2pV2/2e  =  0  (32) 

This  equation  is  well  behaved  and  simple  to  integrate  numerically.  The  boundary  condition 
employed  for  t  is  taken  to  be 

-f-  =  0  (33) 

dy 

While  it  is  realized  that  this  condition  may  not  be  entirely  accurate  at  the  wall  [21],  it  is  used 
at  this  initial  stage  for  simplicity  and  may  be  modified  in  further  work. 

In  regions  removed  from  the  wall,  the  e  equation  remains  well  behaved.  The  k  equation 

d2f 

on  the  other  hand  contains  the  source  term  k  V2/  2.  V2/2,  which  can  be  approximated  as  — ~ 

dy 

has  the  following  properties 


V2/2  >0  y+<  3.4  (34) 

V2/2  <  0  y+  >  3.4 

where  y+  =  3.4  represents  the  point  of  inflection  in  the  /2  function.  In  regions  where  V2/2  is 
negative  or  zero,  the  k  equation  is  well  behaved.  However,  V2/2  large  and  positive  represents  a 
growing  source  term,  which  can  be  numerically  unstable.  However,  since  the  point  y+  =  3.4  is 
very  close  to  the  wall,  and  within  the  laminar  sublayer,  k  can  be  approximated  by  the  relation 


or 


k*  =  constant  .  y+2 


(35) 


k 


=  k 


„+2 


wall 


D  -  exp(-^-)]2 


(36) 


which  from  direct  simulations  [22],  is  generally  known  to  be  valid  up  to  y+  =  10.  Finally,  in 
regions  far  away  from  the  wall,  the  damping  functions  become  unity,  their  derivatives  all  van¬ 
ish,  and  the  original  high-Reynolds-number  equations  are  recovered,  albeit  with  the  new  values 
of  the  constants  advocated  in  [18].  Thus,  in  summary,  the  e  equation  given  in  the  form  (30)  is 
employed  throughout  the  entire  flow-field,  except  at  the  wall,  where  the  form  (32)  is  used.  For 
the  it  equation,  the  form  given  by  equation  (30)  is  employed  from  the  far-field  up  to  y+  =  3.4, 
which  is  within  the  laminar  sublayer.  Below  this  value  of  y+,  equation  (36)  is  employed  with 
the  boundary  condition  for  k  given  by  equation  (31). 


The  multigrid  strategy  previously  described  for  the  high-Reynolds  number  turbulence 
equations  carries  over  in  a  straight-forward  manner.  The  linearization  of  the  source  terms  is 
now  taken  as 


dH 

dw 


2  c  .  V2/2 

-  -S  +  (p  +  —  )-77- 

3  S/t  PJ  2 

ClC^f}J2P+C2^ 


/  2 


-|c,S-2C,f 


(37) 


where  the  production  term  in  the  e  equation  has  been  simplified  by  the  definition  of  (i,  in  equa¬ 
tions  (27),  in  order  to  remove  k  from  the  denominator.  The  damping  functions  are  evaluated 
only  on  the  finest  grid,  and  interpolated  up  to  the  coarser  grids,  thus  affording  a  more  con¬ 
sistent  representation  of  the  equations  on  all  grid  levels.  A  full  multigrid  strategy  is  employed, 
where  grid  sequencing  is  used  to  provide  an  initial  solution  for  the  fine  grid.  In  general,  it  has 
been  found  advantageous  to  use  the  high-Reynolds-number  model  with  wall  functions  on 
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coarse  grids,  and  the  low-Reynolds-number  model  on  the  finest  grid  when  grid  sequencing, 
thus  rapidly  setting  up  appropriate  levels  of  eddy  viscosity  on  the  finest  grid. 

10.  RESULTS 

The  present  implementation  of  the  low-Reynolds-number  turbulence  model  has  been 
employed  to  compute  the  turbulent  boundary  layer  over  a  flat  plate,  the  transonic  flow  over  an 
RAE  2822  airfoil,  and  the  transonic  flow  over  a  two-element  airfoil. 

The  mesh  employed  to  compute  the  flat  plate  boundary  layer  case  is  depicted  in  Figure 

11.  It  contains  24  points  ahead  of  the  plate,  48  points  along  the  plate  in  the  streamwise  direc¬ 
tion,  and  80  points  in  the  direction  normal  to  the  plate.  The  freestream  Mach  number  is  0.3, 
and  the  Reynolds  number  of  the  flow,  based  on  the  length  of  the  plate  is  10  million.  The  first 
point  normal  to  the  plate  is  located  at  a  distance  of  2x  10-6  plate  lengths,  which  lies  in  the 
region  y+<l.  The  resulting  velocity  profiles  are  plotted  in  Figures  12  and  13,  both  in  physical 
coordinates,  and  logarithmic  wall  coordinates,  and  compared  with  the  well  known  l/7th  power 
law  distribution,  and  logarithmic  law  of  the  wall  profile.  The  computed  skin  friction  is  plotted 
in  Figure  14,  versus  the  experimental  data  taken  from  [23].  The  resulting  distributions  of  k 
and  e  are  shown  in  Figures  15  and  16.  The  well  known  peaks  of  jfc  and  e  are  observed,  and  a 
non-zero  value  of  e  at  the  wall  is  obtained.  These  distributions  are  however  slightly  different 
from  those  obtained  previously  with  the  same  model  [18],  and  may  be  attributed  either  to  the 
different  boundary  condition,  or  to  the  near-wall  grid  resolution.  The  overall  flow  quantities  are 
nevertheless  well  predicted,  as  shown  in  Figures  12  through  14. 

The  transonic  flow  case  over  the  RAE  2822  airfoil  presented  in  the  previous  section  has 
been  recomputed  with  the  low-Reynolds-number  turbulence  model  (Mach  =  0.729,  Incidence  = 
2.31  degrees.  Re  =  6.5  million).  The  mesh  employed  is  similar  to  that  shown  in  Figure  1, 
except  that  the  normal  spacing  at  the  wall  is  now  reduced  to  lx  10-5  chords,  which  results  in 
the  first  mesh  point  off  the  wall  in  the  region  1  <  y+  <  3  over  the  entire  surface  of  the  airfoil. 
The  computed  Mach  contours  and  eddy  viscosity  contours  are  similar  to  those  depicted  in  Fig¬ 
ures  2  and  3,  except  in  the  near-wall  regions,  where  both  quantities  vanish  rapidly.  The  com¬ 
puted  surface  pressure  and  skin-friction  distributions  are  compared  with  experimental  data  in 
Figures  17  and  18.  The  computed  lift  is  slightly  lower  than  that  predicted  with  the  wall  func¬ 
tions  and  that  previously  obtained  using  an  algebraic  model  [7].  At  present,  it  is  not  clear 
whether  this  is  due  to  the  actual  model  formulation,  or  is  associated  with  the  present  imple¬ 
mentation  (artificial  dissipation,  grid  resolution).  However,  the  differences  are  rather  small  and 
the  skin  friction  appears  to  be  well  predicted.  The  convergence  of  the  density  equation  and  the 
two  turbulence  equations  is  depicted  in  Figure  19,  where  the  residuals  are  plotted  versus  the 
number  of  multigrid  cycles  on  the  finest  grid.  The  flowfield  and  turbulence  equations  are  all 
initialized  with  uniform  freestream  values,  and  25  cycles  were  performed  on  the  previous 
coarser  grid  using  wall  functions,  prior  to  initializing  the  solution  procedure  on  the  finest  grid. 
Initializing  the  calculation  with  freestream  values  for  all  equations  on  the  finest  grid  has  also 
been  employed  with  little  degradation  in  convergence.  From  Figure  19,  all  equations  are  seen 
to  converge  at  approximately  the  same  rate,  resulting  in  a  residual  reduction  of  4  to  5  orders  of 
magnitude  over  300  multigrid  cycles. 

The  final  case  involves  the  transonic  flow  over  a  two-element  airfoil.  This  case  illustrates 
the  ease  with  which  complex  geometries  and  flows  with  multiple  viscous  layers  may  be  han¬ 
dled  by  the  present  methodology.  The  mesh  employed  is  depicted  in  Figure  20.  It  contains  a 
total  of  28,871  vertices,  with  a  normal  spacing  of  2x1  Or5  chords  off  the  wall  for  each  airfoil 
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element.  The  freestream  Mach  number  is  0.5,  the  incidence  is  7.5  degrees,  and  the  Reynolds 
number  is  4.5  million.  The  computed  Mach  contours  and  eddy  viscosity  contours  are  depicted 
in  Figures  21  and  22.  At  these  conditions,  the  flow  is  supercritical  and  a  shock  is  formed  on 
the  upper  surface  of  the  slat  A  small  region  of  separated  flow  occurs  behind  the  shock,  as  can 
be  seen  from  the  skin  friction  plot  of  Figure  23.  This  region  of  separation  has  previously  been 
reported  in  prior  calculations  using  an  algebraic  turbulence  model  [7].  The  computed  surface 
pressure  distribution  is  seen  to  compare  favorably  with  experimental  wind-tunnel  data  [24],  in 
Figure  24.  The  convergence  rate  for  this  case  is  depicted  in  Figure  25,  where  the  residuals  of 
the  density  equation  and  the  two  turbulence  equations  are  reduced  by  approximately  3  to  4  ord¬ 
ers  of  magnitude  over  300  cycles  on  the  finest  grid. 

11.  CONCLUSION 

A  multigrid  strategy  for  solving  the  steady-state  high  and  low-Reynolds  number  k  -  e  tur¬ 
bulence  equations  has  been  formulated  and  implemented  on  unstructured  meshes.  A  variety  of 
aerodynamic  flows  have  been  computed,  consistently  demonstrating  similar  convergence  rates 
for  the  turbulence  and  flow  equations.  Initialization  of  all  flow  and  turbulence  quantities  may 
be  performed  using  uniform,  freestream  values.  At  present,  the  evaluation  of  the  turbulence 
terms  requires  a  significant  fraction  of  the  overall  time  within  a  single  time-step.  For  example, 
the  RAE  2822  supercritical  airfoil  flow  case  with  the  low-Reynolds-number  turbulence  model 
requires  roughly  2.5  seconds  per  multigrid  cycle  on  a  single  processor  of  the  CRAY-YMP 
supercomputer,  which  is  almost  75%  higher  than  that  required  by  the  algebraic  model  reported 
previously  [7].  However,  it  is  estimated  that  this  can  be  substantially  reduced  by  assembling 
the  turbulence  and  flow  residuals  simultaneously  within  a  single  loop.  Given  the  demonstrated 
convergence  rates,  the  two-equation  turbulence  model  should  be  competitive  in  terms  of  com¬ 
puter  resources  with  algebraic  models,  while  providing  much  greater  flexibility  in  dealing  with 
complex  geometries  and  flow- fields. 

Future  work  should  involve  a  more  thorough  investigation  of  the  various  two-equation 
turbulence  models  and  their  ability  in  predicting  complex  aerodynamic  flows,  including  flows 
with  massive  separation. 
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Figure  1 

Unstructured  Mesh  Employed  for  Computing  Flow  Over  an  RAE  2822  Airfoil 
(Number  of  Vertices  =  12,823) 


Figure  2 

Computed  Mach  Contours  for  Turbulent  Flow  over  RAE  2822  Airfoil 
Using  Wall  Functions  (Mach  =  0.729,  Re  =  6.5  million.  Incidence  =  2.31  degrees) 


Figure  3 

Computed  Eddy  Viscosity  Contours  for  Turbulent  Flow  over  RAE  2822  Airfoil 
Using  Wall  Functions  (  Mach  =  0.729,  Re  =  6.5  million,  Incidence  =  2.31  degrees ) 
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Figure  5 

Convergence  Rale  of  Density  Equation  and  K  equation  Using  Wall  Functions 
versus  the  Number  of  Multigrid  Cycles  for  Flow  over  an  RAE  2822  Airfoil 
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Figure  6 

Unstructured  Mesh  Employed  for  Computing  Flow  Over  a  Four-Element  Airfoil 
(Number  of  Vertices  =  51,100) 


Figure  7 

Computed  Mach  Contours  Using  Wall  Functions  for  Flow  ova  a  Four-Element  Airfoil 
(Mach  ■  0.2,  Re  *  2.83  million,  Incidence  =  8.18  degrees) 
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Figure  9 

Comparison  of  Computed  Surface  Pressure  using  Wall  Functions  with  Experimental 
Wind-Tunnel  Data  for  Flow  over  a  Four-Element  Airfoil 
(Mach  =  0.2,  Re  =  2.83  million,  Incidence  =  8.18  degrees) 
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Figure  11 

Triangular  Mesh  Employed  for  Flat  Plate  Boundary  Layer  Calculation 
(Number  of  Vertices  =  5913,  10:1  Magnification  in  Y-direction) 
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Figure  12 

Computed  Velocity  Profile  in  Physical  Coordinates  Versus  the  l/7th  Power  Law  Profile 

(Mach  =  0.3,  Re,  =  5.3  million) 
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Figure  13 

Computed  Velocity  Profile  in  Logarithmic  Coordinates  Versus  Logarithmic  Law  of  the  Wall 

(Mach  =  0.3,  Re,  =  5.3  million) 


0.00  1.00  ZOO  3.00  4.00  5.00  6.00  7.00  g.00  9.00  10.00 


Figure  14 

Computed  Skin  Friction  Distribution  for  Flat  Plate  Boundary  Layer 
Versus  Experimental  Data  from  [23] 
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Figure  16 

Computed  Near  Wall  Distribution  of  e*  for  Flat  Plate  Boundary  Layer 
(Mach  =  0.3,  Re,  =  3.3  million) 
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Figure  17 

Computed  Surface  Pressure  Distribution  Using  Low-Reynolds  Number  Modification 
for  Turbulence  Equations  Versus  Experimental  Data  for  Flow  past  RAE  2822  Airfoil 
(Mach  ■  0.729,  Re  =  6.S  million,  Incidence  *  2.31  degrees) 
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Figure  18 

Computed  Skin-Friction  Distribution  Using  Low-Reynolds  Number  Modification 
for  Turbulence  Equations  Versus  Experimental  Data  for  Flow  past  RAE  2822  Airfoil 
(Mach  =  0.729,  Re  =  6.5  million,  Incidence  =  2.31  degrees) 
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Figure  19 

Convergence  Rate  of  the  Density  Equation  and  the  Two  Turbulence  Equations 
Modified  for  Low-Reynolds  Number  Effects  Versus  the  Number  of  Multigrid  Cycles  for  Flow  Past  RAE  2822  Airfoil 
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Flgnre20 

Global  View  of  Come  Unstructured  Mesh  aid  Close-Up  View  of  Fine 
Unstructured  Merit  Employed  for  Computing  Flow  Put  a  Two-Element  Airfoil 
(Come  Mesh  Points  •  7272,  Fine  Mesh  Points  ■  28871) 


Figure  21 

Computed  Mach  Contours  Using  Low-Reynolds  Number  Modification  for  Turbulence 
Equations  for  Supercritical  Flow  over  a  Two-Element  Airfoil 
(Mach  *  0.3,  Re  ■  4 3  million.  Incidence  ■  7  J  degrees) 


Figure  22 

Computed  Eddy  Viscosity  Contours  Using  Low-Reynolds  Number  Modification 
for  Turbulence  Equations  for  Supercritical  Flow  over  a  Two-Element  Airfoil 
(Mach  =  0.5,  Re  =  4.5  million.  Incidence  -  7.5  degrees) 
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Figure  23 

Computed  Skin-Friction  Distribution  on  Slat  Showing  Region  of  Separated 
Flow  Behind  Upper  Surface  Shock  (Mach  =  0.5,  Re  =  4.5  million.  Incidence  =  7.5  degrees) 


Figure  24 

Computed  Surface  Pressure  Distribution  Using  Low-Reynolds  Number  Modification 
for  Turbulence  Equations  for  Supercritical  Flow  over  a  Two-Element  Airfoil 
Versus  Experimental  Wind  Tunnel  Data 
(Mach  =  0.5,  Re  =  4.5  million,  Incidence  =  7.5  degrees) 
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Figure  25 

Multigrid  Convergence  Rate  of  the  Density  Equation  and  the  Two  Turbulence 
Equations  Using  Low-Reynolds  Number  Modifications  for  Flow  Over 
Two-Element  Airfoil  (Mach  =  0.5,  Re  =  4.5  million.  Incidence  =  7.5  degrees) 
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